/* Table 1. Descriptives */ 
proc means data=jcms2023 mean std var range min max kurtosis skewness n;
 var count staff_influence membe_influence access_strategies q21_01n staffsize age agelog breath;

proc freq data=access; 
 tables membe_influence*revolvingdoor/chisq;

proc freq data=access; 
 tables count newtype revolvingdoor/missing;

/* imputation of missing values */
proc mi data=access;
 em out=test;
 var count staff_influence membe_influence staffsize access_strategies agelog breath;

/* some descriptive results */
proc corr;
 var count staff_influence membe_influence access_strategies staffsize age agelog breath;

proc means mean std maxdec=2;
 var staff_influence;
 class revolvingdoor;

/* categorical variable for figures, plus creating smaller dataset focused on analyses */
data test (keep=cigid staff_influence Staffinfluence revolvingdoor revolving newtype count membe_influence staffsize q21_01n newtype access_strategies breath agelog breath age);
 set test;

if staff_influence le 0.68                             then Staffinfluence = 'Low';
if staff_influence gt 0.68 and staff_influence le 0.80 then Staffinfluence = 'Medium';
if staff_influence gt 0.80                             then Staffinfluence = 'High';

attrib revolving length=$25;;
if revolvingdoor eq '1' then revolving = 'Revolving';
if revolvingdoor eq '0' then revolving = 'Not revolving';

/* zero-inflated poisson model */
proc genmod data=test plots=all;
 class revolvingdoor(ref='0') newtype(ref='1');
 model count=revolvingdoor|staff_influence revolvingdoor|membe_influence staffsize newtype access_strategies breath agelog breath/dist=zip type3 wald;
 zeromodel revolvingdoor|staff_influence revolvingdoor|membe_influence staffsize newtype access_strategies breath agelog breath; 
 output out=zip predicted=pred;
 store out=zip1;
 ods output Modelfit=fit;

/* scoring confidence intervals */
proc plm source=zip1;
 score data=test out=plmout pred stderr lclm uclm/ilink;

/* figure 2 */
proc means data=plmout;
 var Predicted LCLM UCLM;
 class Staffinfluence Revolvingdoor;
*where count gt 0;

/* this prints/exports the data which is then transported to excel */
proc export data=plmout outfile="Z:\3_5_population and CIG-survey project\2_publications\10_JCMS 2022\drafts\submission JCMS\minor revisions\predicted.xlsx" 
            dbms=excel label replace;
     sheet="merge"; 
     newfile=yes;

proc plm restore=zip noinfo;
   effectplot interaction(x=staff_influence sliceby=revolvingdoor) / clm;
/* effectplot interaction(x=Origin sliceby=Type) / clm connect; */ /* or connect the means */

/* figure 3 */
data zip;
 set zip;

logpredicted=log(pred+1);
label logpredicted    ='Log predicted number of meetings';
label staff_influence ='Staff influence';

proc sgplot data=zip noautolegend;
  reg  x=staff_influence y=logpredicted/degree=1 group=revolving clm='95% Confidence Limits';  
  yaxis label='Log predicted number of meetings'; 
  Xaxis label='Staff influence'; 
  keylegend / location=outside position=bottom across=1;
* where logpredicted gt 0;
run;

